
%[0.0000000000000000000002]
%[0.00000000000002]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Paramater
%% Keep the same value in Simulink

global Real_delay_ns;
global Real_delay;
%Real_delay=0.1e-9;  %units of s

Real_delay=0.1e-9;  %units of s

Input_A=1e-3; %units of W=J/s

Samplying_time=0.00000000005; %units of s
Samplying_N=3e-6/Samplying_time

Real_delay_ns=Real_delay*1.0e9; %units of ns
%%%%%%%%%%%%%%%%%%%%%%%%%%%

for C_i=1:30
%for C_i=11:11

    Input_A=10^(-7+(C_i-1)*4/30);
    Matlab_input_Pin=Input_A;


    myVar = Simulink.Parameter;
    myVar.Value = Input_A;           %
    myVar.CoderInfo.StorageClass = 'ExportedGlobal';
    out= sim('Hjh_simulation_withPD_low_power','StartTime', '0', 'StopTime', '0.0003'); %100 measurement
   %out= sim('Hjh_simulation_withPD_low_power','StartTime', '0', 'StopTime', '0.0015'); %500 measurement
   %out= sim('Hjh_simulation_withPD_low_power','StartTime', '0', 'StopTime', '0.00006');  %20 measurement
   

    %%%%%% calculate CFI information
    I_normnal=0.0;
    for i=1:60000
        tt=Samplying_time*i;
        I_normnal=I_normnal+Gaussian_pulse_delay(tt,Real_delay);
    end
    % I_normnal

    %%%%%% calculate CFI information
    I=0.0;
    for i=1:60000
        tt=Samplying_time*i;
        I=I+(1/I_normnal)*Gaussian_pulse_delay(tt,Real_delay)*( ( log10((1/I_normnal)*Gaussian_pulse_delay(tt,Real_delay))-log10((1/I_normnal)*Gaussian_pulse_delay(tt,0.0)) )/Real_delay )^2;
    end

    % Physical constants
    h = 6.626e-34;      % Planck's constant (J·s)
    c = 3e8;            % Speed of light (m/s)
    lambda = 633e-9;    % Wavelength (m)

    % Energy per photon
    E = h * c / lambda;  % Unit: joules (J)

    % Convert optical power (W) to photon count per second
    photons = max(Input_A, 0)*3e-6*1*20 / E % 3e-6 s is detection time; 1 is postselection;
    CFI= photons*I;
    CFI_Std= 1e9*1/sqrt(abs(CFI)) %units of ns

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Data process
    x1=out.APD.signals(1).values;
    TT=60000;
    [peaks, locs1] = findPeriodicPeaks(x1, TT);
    %disp('Locs1： ');   disp(locs1);
    %disp('Peaks：');   disp(peaks);

    x2=out.APD.signals(2).values;
    [peaks, locs2] = findPeriodicPeaks(x2, TT);
    %disp('Locs2： ');   disp(locs2);

    Delta_t_WVA_M1=abs( (std(locs1-locs2))  *Real_delay_ns/(mean(locs1-locs2)) )%units of ns

    samples_per_period = TT;%

    y=out.ACI.signals(1).values;
    num_periods = floor(length(y) / samples_per_period);
    period_max = zeros(1, num_periods);

    for i = 1:num_periods
        start_idx = (i-1)*samples_per_period + 1;
        end_idx = i*samples_per_period;
        period_max(i) = max(y(start_idx:end_idx));
    end
    peaks_ACI0=period_max;


    y=out.ACI.signals(2).values;
    num_periods = floor(length(y) / samples_per_period);
    period_max = zeros(1, num_periods);

    for i = 1:num_periods
        start_idx = (i-1)*samples_per_period + 1;
        end_idx = i*samples_per_period;
        period_max(i) = max(y(start_idx:end_idx));
    end
    peaks_ACI=period_max;


    peaks_ACI0 = peaks_ACI0(2:end);
    peaks_ACI = peaks_ACI(2:end);

    mean(peaks_ACI);
    mean(peaks_ACI0);

    Delta_t_AWVA=abs( (std(peaks_ACI-peaks_ACI0))  *Real_delay_ns/(mean(peaks_ACI-peaks_ACI0)) )%units of ns



    fid = fopen('Precision_dependence_of_photons_low_power.txt', 'a');
    %fprintf(fid,'# Input_A  CFI WVA AWVA\n');
    fprintf(fid,'%e %e %e %e\n',...
        Input_A*1000,...  % mW
        CFI_Std,...
        Delta_t_WVA_M1,...
        Delta_t_AWVA...
        );
    fclose(fid);

end
